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Abstract. A discussion of methods for reducing the noise variance of flavor singlet 
quantities ("disconnected diagrams") in lattice QCD is given. After an introduc- 
tion, the possible advantage of partitioning the Wilson fermion matrix into disjoint 
spaces is discussed and a numerical comparison of the variance for three possible 
partitioning schemes is carried out. The measurement efficiency of lattice operators 
is examined and shown to be strongly influenced by the Dirac and color partition- 
ing choices. Next, the numerical effects of an automated subtraction algorithm on 
the noise variance of various disconnected loop matrix elements are examined. It is 
found that there is a dramatic reduction in the variance of the Wilson point-split 
electromagnetic currents and that this reduction persists at small quark mass. 



1 Introduction 

1.1 Motivations 

The calculation of flavor singlet quantities, also referred to as disconnected 
diagrams because the fermion lines are disjoint, is one of the greatest technical 
challenges left in lattice QCD. Disconnected contributions are present in a 
wide variety of quantities in strong interaction physics including all baryon 
form factors, axial operators involving quark spin context (Ellis- Jaffe sum 
rule), hadronic coupling constants and polarizabilities, the p-N sigma term, 
and propagation functions for various flavor singlet mesons. Such quantities 
are also present in deep inelastic structure functions measured on the lattice 
using the operator product expansion, but are not included because of their 
difficulty and large Monte Carlo error bars. These types of diagrams are 
difficult to evaluate because exact extractions require many matrix inversions 
to measure all the background fermion degrees of freedom (including space- 
time). Disconnected quark contributions are instead isolated stochastically 
by a process of applying "noises" to the fermion matrix to project out the 
desired operator contribution. For an overview of selected aspects of flavor 
singlet calculations in lattice QCD, see Ref. [Q. 

1.2 Mathematical Background 

Noise methods are based upon projection of the signal using random noise 
vectors as input. That is, given 
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Mx = rj, 

where M is the quark matrix, x is the solution vector and 77 is the noise 
vector, with 

< ry^ >= 0, < rj^rij >= (5y, 

where one is averaging over the noise vectors, any inverse matrix element, 
M,7^, can then be obtained from 

< ^ox, < VjVk >= M-^. 

k 

We shall consider two techniques for reducing the noise variance in lattice 
QCD simulations: partitioning[Q and subtraction methods]^. Partitioning 
the noise appropriately, which means most generally a zeroing out of some 
pattern of noise vector elements, but which will specifically be implemented 
here by using separate noise source vectors in Dirac and color spaces, can 
lead to significant reductions in the variance. Subtraction methods, which 
involve forming new matrix operators which have a smaller variance but the 
same expectation value (i.e. are "unbiased"), can also be of great help. The 
key is using a perturbative expansion of the quark matrix as the subtraction 
matrices which, however, are not unbiased in general and require a separate 
calculation, either analytical or numerical, to remove the bias. We will see 
that the various lattice operators have dramatically different behaviors under 
identical partitioning or subtraction treatments. Both of these methods, par- 
titioning and subtraction, will be treated numerically but the hope is that the 
numerical results will eventually be "explained" by some simple rules based 
on the structure of the Wilson matrix. 



2 Noise Theory 

2.1 Variance Evaluations 

Let us review the basics of matrix inversion using noise theory. The theoret- 
ical expressions for the expectation value and variance (V) of matrices with 
various noises are given in Rcf. Q| . One has that 



1 ^ 

Xynn = Y^^rinillfnl- (1) 



1=1 



(m, n = 1, . . . ,N] I — 1, . . . , L.) We have 
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and the expectation value, 

By definition the variance is given by 

V[Tr{QX}] =< I J2 ImnXnm - Tr{Q}\'' > . (4) 

The variance may be evaluated as, 

V [Tr{QX}] =J2(< l^""'!" > (5) 

+ qrnnqnm < {Xmnf >) + < - Ip > knnl^ ■ 

n 

2.2 Real Noises 

For a general real noise, 

< l^mnP >= J, (6) 
< (Xmn)' >= (7) 

for TO ^ n so that 

F[rr{QX,.eai}] = ;^ E d?""'!' + 9™"'?mJ (8) 

n 

The case of real Z{2) has Eqs.(|) and (|) holding for m =^ n, but also 

< |X„„ - 1|2 >=0. (9) 

This shows that 

V[Tr{QXzi2)}] < V[Tr{QX,,,i}]. (10) 
Thus, Z{2) noise has the lowest variance of any real noise. 

2.3 General Z(N) Noise 

For general Z[N) [N > 3) noise we have a different situation. One has that 

< 1 (11) 

< {X^nf >- 0, (12) 
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for m ^ n, and again 

< |X„„-1|2 >=0. (13) 

Thus 

V[Tr{QXz(N)}] = \Y. (14) 

and the variance relationship of Z{2) and Z{N) is not fixed for a general 
matrix Q. The reason for the difference in Eqs.(Q) and ([l^) is that the square 
of an equally weighted distribution of Z{2) elements is not itself uniformly 
distributed (always 1), whereas the square of a uniformly weighted Z{N) 
distribution for > 3 is also uniformly distributed. However, if the phases 
of qmn and 9*„ are uncorrelated, then V\Tr{QXz{2)}] « y[Tr{QXz[N)}], 
{N > 3) which, we will see, is apparently the case for the operators studied 
here. 



3 Partitioning the Problem 
3.1 Basic Idea 

By "partitioning" I mean replacing the single noise vector problem, 

Mx = rj, (15) 

which yields a complete output column, M^rjk, from a single input noise 
vector, 77, with a problem 

MxP = r]P,p^l,...P, (16) 

where the 77^' have many zeros corresponding to some partitioning scheme. 
In this latter case it takes P inverses to produce a complete measurement or 
sampling of a column of M^^. 

For the unpartitioncd problem (for Z{N), N > 3, say) 

V[Tr{QX}] = 1 ^ |g„„p, (17) 

= ^NiN-l)<\q\^>, 

where I have defined the average absolute squared off diagonal matrix ele- 
ment, < |gp >. For the partitioned problem the total variance includes a 
sum on p, 

j2V[Tr{QX,}] = ^J2 E I'?".".!'' (18) 
^l7V(^-l)<|gpp>. 



Noise Methods 5 



In order for this method to pay off in terms of computer time, one needs 
that 

J2 V[Tr{QXp}] < -V[Tr{QX}], (19) 
p=i 

< |9P > . (20) 

The goal of partitioning is to avoid some of the large off-diagonal matrix 
elements so that in spite of doing P times as many inverses, a smaller variance 
is produced for the same amount of computer time. The spaces partitioned 
can be space-time, color or Dirac or some combination. I have found that 
partitioning in Dirac and color spaces can strongly affect the results. 

3.2 Simulation Description 

I consider all local operators, '4)(x)r'il){x), as well as point-split versions of the 
vector and axial vector operators. This means 16 local operators and 8 point- 
split ones, making a total of 24, which are listed below. For each operator 
there are both real and imaginary parts, but in each case one may show via the 
quark propagator identity S = 756'''' 75, that only the real or the imaginary 
part of each local or nonlocal operator is nonzero on a givc;n configuration for 
each space-time point. However, this identity is not respected exactly by noise 
methods, so the cancellations are actually only approximate configuration by 
configuration. However, the knowledge that one part is purely noise allows 
one to simply drop that part in the calculations, thus reducing the variance 
without biasing the answer. 
The operators I consider are: 

• Scalar: Re[^l){x)il'{x)\ 

• Local Vector: /TO[^/i(x)7^V(2;)] 

• Point-Split Vector: 

Klm[i){x + a^)(l +j^)Ul{x)ip{x) - V'(a;)(l - 7^)f/^(a;)V'(a; -|- a^)] 

• Pseudoscalar: Re[^p{x)j5'tp{x)] 

• Local Axial: Re['tp{x)^5jij,tlj{x)] 

• Point-Split Axial: 

KRe[4>{x + a^_)757^f/^(x)V'(a;) + Vi(a;)757^Z7^(a;)V'(a; -|- a^)] 

• Tensor: Im[ip{x)a fj,,y'ip{x)] 

I actually consider the zero momentum version of these operators, summed 
over both space and time. 

The sample noise variance in M quantities Xi is given by the standard 
expression: 

Vnoise = - (21) 

i=l,M 
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What I concentrate on here are the relative variances between the different 
methods. Since the squared noise error in a single configuration is given by 

the ratio of variances gives a direct measure of the multiplicative ratio of 
noises, and thus the computer time, necessary to achieve the same noise error. 
However, the variance itself does not take into account the extra P inverses 
done when problem is partitioned. In order to measure the relative efficiency 
of different partitionings, I form what I call pseudo-efficiencies ratios ("PE"), 
which are defined by 

methodl _ P.,nethodl{Vno.seT''''°''^ . . 

Sneth0d2^ " P™et/.orfl(Ko»se)™^*'*°'*2' ^ ) 

where Pmethod are the number of partitions required by the method. I refer 
to these ratios as "pseudo" efficiencies since I do a fixed number of iterations 
for all of the operators I consider. It could very well be that different methods 
will require significantly different numbers of iterations of conjugate-gradient 
or minimum residual for the same level of accuracy. One desires to find the 
lowest PE ratio for a given operator. |^ 

I display PE ratio results for Wilson fermions in a 16^ x 24, (3 = 6.0 lattice 
with K = 0.148 in Table 1, which follows on the next page. (Part of this Table 
also appeared in Ref. (P.) I will examine 3 partitionings: 

• Z{2) unpartitioned {''P ^ 1 Z(2)"); 

• Z{2) Dirac partitioned ("P = 4 Z(2)"); 

• Z{2) Dirac and color partitioned ("P = 12 ^(2)"). 

My Dirac gamma matrix representation is: 

= {I ) '^'^ = (o -i) = " ^0 ■ ^^^^ 

In Table 1 I list the relative PEs of the two partitioned methods rel- 
ative to the unpartitioned case. Referring to the above list of the real or 
imaginary parts of operators my notation here is as follows: "Scalar" stands 
for the operator "Local Vector 1" for example stands for the operator 
'i/'7i'0,"P-S Vector 1" stands for the 1 component of the point split vector 
current, "pseudoscalar" stands for Vj75'(/;, "Local Axial 1" stands for the oper- 
ator ■!/;757i^,"P-S Axial 1" stands for the point split axial 1 component, and 
"Tensor 41" stands for example for the operator tpaui/j. 

^ Of course an evaluation via TV partitioning on an x A'^ matrix yields a PE 
numerator factor of zero relative to other methods since the variance is exactly 
zero in this case. This is of course prohibitively expensive, but it suggests that 
more efficient partitionings are possible for large computer budgets. Thanks to 
M. Peardon for bringing this point out. 
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There are extremely large variations in the behaviors of the operators 
listed in Table 1 under identical partitionings. Of the partitionings considered 
it is most efficient to calculate scalar and vector operators with an unparti- 
tioned simulation. On the other hand, it is far more efficient to calculate the 
pseudoscalar in a Dirac and color partitioned manner. Notice the entries for 
the 1,2 components of the axial current (both local and point split) do not 
behave like the 3,4 components under pure Dirac partitioning, but they do 
when Dirac and color partitionings are combined. Four of the tensor opera- 
tors respond best to a pure Dirac partitioning, while the other two prefer a 
partitioning in Dirac and color spaces combined. The ratio of the largest to 
the smallest entry in the right hand column is about 800! 

As pointed out in Section 2, the variance of Z{2) and Z{N) {N > 3) noises 
are in general different. For this reason I also investigated partitioning using 
Z(4) as well as volume (gauge variant) noises, but there do not seem to be 
large factors to be gained relative to the Z{2) case. The SESAM collaboration 
also has seen the efficacy of partitioning (in Dirac space) for axial operators B . 



Table 1. The pseudoefficiency (PE) ratios associated with the methods in- 
dicated. 

Operator PE(gif|§) PEl ^p-^"'^' ) 



Scalar 


2.83 


± 


0.47 


10.9 


± 


2.5 


Local Vector 1 


2.38 


± 


0.65 


8.71 


± 


1.8 


Local Vector 2 


2.50 


± 


0.53 


12.1 


± 


2.9 


Local Vector 3 


3.60 


± 


1.00 


11.4 


± 


2.4 


Local Vector 4 


3.41 


± 


0.60 


16.3 


± 


3.4 


P-S Vector 1 


2.63 


± 


0.56 


9.94 


± 


2.2 


P-S Vector 2 


2.27 


± 


0.44 


11.0 


± 


2.3 


P-S Vector 3 


3.52 


± 


0.74 


11.4 


± 


1.5 


P-S Vector 4 


3.87 


± 


0.49 


15.5 


± 


4.2 


Pseudoscalar 


0.698 


± 


0.15 


0.0201 


± 


0.0043 


Local Axial 1 


0.114 


± 


0.021 


0.144 


± 


0.029 


Local Axial 2 


0.126 


± 


0.020 


0.146 


± 


0.037 


Local Axial 3 


1.13 


± 


0.19 


0.162 


± 


0.038 


Local Axial 4 


2.26 


± 


0.24 


0.187 


± 


0.035 


P-S Axial 1 


0.167 


± 


0.040 


0.151 


± 


0.018 


P-S Axial 2 


0.110 


± 


0.032 


0.114 


± 


0.028 


P-S Axial 3 


1.67 


± 


0.35 


0.186 


± 


0.049 


P-S Axial 4 


1.85 


± 


0.21 


0.238 


± 


0.036 


Tensor 41 


1.07 


± 


0.30 


0.295 


± 


0.049 


Tensor 42 


0.345 


± 


0.076 


0.0889 


± 


0.011 


Tensor 43 


1.32 


± 


0.43 


0.398 


± 


0.13 


Tensor 12 


1.12 


± 


0.25 


0.376 


± 


0.066 


Tensor 13 


0.116 


± 


0.024 


0.363 


± 


0.053 


Tensor 23 


0.0314 


± 


0.0058 


0.0751 


± 


0.016 
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4 Perturbative Noise Subtraction 

4.1 Description of Algorithm 

Consider Q such that 

< Tr{QX} >= 0. (25) 

One can then form 

<Tr{{Q-Q)X}>=<Tr{QX}> . (26) 

However, 

V[Tr{{Q - Q)X}] ^ V[Tr{QX}]. (27) 

As we have seen for Z{N) {N > 2), the variance comes exclusively from off 
diagonal entries. So, the trick is to try to find matrices Q which are traceless 
(so they do not affect the expectation value) but which mimic the off-diagonal 
part of Q as much as possible to reduce the variance. 

The natural choice is simply to choose as Q the perturbative expansion of 
the quark matrix. Formally, one has (/, J = {x, a, a}) 

iM-')ij = T (28) 

dij - kPij 

where 

Pi' = + (1 - If^Pli^ - a^^)S..y+aJ. (29) 

A' 

Expanding this in k gives the perturbative (or hopping parameter) ex- 
pansion. 



M-^ = I +kP + K^P^ + K^P^ + ■ ■■ . (30) 

One constructs < rjj{M^^)ikr]k > and subtracts it from < rjjM^rjk >, 
where 77 is the noise vector. Constructing < r]j{Alp^)ikr]k > is an iterative 
process and is easy to code and extend to higher powers on the computer. I 
will iterate up to 10th order in k. 

One can insert coefficients in front of the various terms and vary them 
to find the minimum in the variance, but such coefficients are seen to take 
on their perturbative value, at least for high order expansions[^. However, 
see also Ref.Q] where in low orders there is apparently an advantage to this 
procedure. Interestingly, significant subtraction improvements occur in some 
operators even in 0th order (point split vectors and two tensor operators.) 

For a given operator, O, the matrix OMp^ encountered in the context of 
< ipOip >= —Tr{OM^^) is not traceless. In other words, one must re-add 
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the perturbative trace, subtracted earher, to get the full, unbiased answer. 
How does one calculate the perturbative part? The exact way is of course to 
explicitly construct all the gauge invariant paths (up to a given k order) for 
a given operator. Another approach is to subject the perturbative contribu- 
tion to a separate Monte Carlo estimation. This is the approach taken here. 
Local operators require perturbative corrections starting at 4th order (ex- 
cept a trivial correction for ipil) at zeroth order) and point split ones require 
corrections starting instead at 3rd order. Because one is removing the bias 
(perturbative trace) by a statistical method, I refer to this as a "statistically 
unbiased" method. Some efficiency considerations in carrying out this proce- 
dure will be discussed in Section 5. Other versions of subtraction methods in 
the context of lattice evaluations of disconnected diagrams may be found in 
Refs.§ and §. 

4.2 Numerical Results 

I am carrying out this numerical investigation in an unpartitioned sense 
(P = 1). The operators which respond best to this partitioning, as discussed 
previously, are the scalar and local and point-split vector currents, and atten- 
tion will be limited here to these cases. The effect of combining the partition 
and subtraction methods has not yet been investigated. I show the ratio of 
unsubtracted variance divided by subtracted variance, Vunsub/Vgub in Figs. 
1 and 2. Factors larger than one give the multiplicative gain in computer time 
one is achieving. The lattices are again Wilson 16'^ x 24, j3 — 6.0. 

Notice the approximate linear rise in the variance ratio as a function 
of subtraction order for the point-split vector charge density at both k = 
0.148 and 0.152, Figs. 1 and 2 respectively. Also notice that even at 5 = 
(subtracting the Kronecker delta) there is a reduction in the variance. The 
slope of the subtraction graph at k = 0.148 is about 3.5; the slope at k = 0.152 
is reduced to a little under 3.0. Although I do not show the results here, the 
same linear behavior is evident in xptp and the local vector operators although 
their slopes are considerably smaller. 

My final results are summarized in Fig. 3, which gives the reduction in 
the variance in the scalar and vector operators after a 10th order subtraction 
has been made at k = 0.148. It is not known why the point split vector 
current responds the best to subtraction. The 10th order point-split vector, 
local vector, and scalar variance ratios change from ~ 35, ~ 12, and ^ 10 at 
K, = 0.148, to ~ 25 ~ 10, and ~ 5 at k = 0.152, respectively. These are all 
zero momentum operators. Although the results are not shown, I have found 
essentially identical results to the above for momentum transformed data, 
necessary for disconnected form factors. Perturbative subtraction methods 
will thus be extremely useful in lattice evaluations of nucleon strangeness 
form factors using the point split (conserved) form of the vector current. 
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Fig. 1. Effect of the level of perturbative subtraction, up to tenth order in k, 
on the ratio of unsubtracted divided by subtracted noise variance for the zero 
momentum point-split (conserved) charge density operator, J4, at k = 0.148. 



Fig. 2. Same as Fig. 1 but for k = 0.152. 
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Fig. 3. Graphical presentation of the unsubtracted variance divided by sub- 
tracted variance of 9 different lattice operators after tenth order subtrac- 
tion at K = 0.148 for Wilson fermions. (Abbreviations used for operators: 
S=scalar; L VEC = local vector; P-S VEC = point split vector). 

5 Efficiency Considerations 
5.1 Fixed Noise Case 

Let me close with some simple observations regarding statistical errors in 
flavor singlet Monte Carlo simulations. There are two sources of variance in 
such simulations: gauge configuration and noise. Given N configurations and 
M noises per configuration, the final error bar on a given operator is given 



where Vg^age and Vnoise are the gauge configuration and noise variances. For 
fixed NM (total number of noises), it is clear that Eq.(^l|) is minimized by 
taking M = I. Thus, in this situation it is best to use a single noise per 
configuration. This simple result is modified by real world considerations of 
overheads. For example, if one assumes that there is an overhead associated 
with generating configurations and fixes instead the total amount of computer 
time, 



by 




T = NM + GnN, 



(32) 



12 Wilcox 



where Gjv is the appropriately scaled configuration generation time overhead, 
then one finds instead that 

M = |^V^, (33) 

>^ gauge 

is the best choice. Note that the ratio Snoise/ S gauge can have a wide range 
of values for various operators, and one is no longer guaranteed that M = 1 
is optimal. 



5.2 Fixed Configuration Case 

Another common real world situation is where N, the number of configura- 
tions, is fixed. In the context of the perturbative subtraction algorithm, one 
should now maximize the number of effective noises for a given computer 
budget. The effective number of noises is 

Meff ^ M{X + SAs), (34) 

where M^jf replaces M in Eq.(|3l|). (M retains its meaning as the actual 
number of gauge field noises.) S is the subtraction order, ranging from to 
10 in Figs. 1 and 2, and As is the slope. X is the factor one obtains from this 
method at 5 = 0, without extra overhead. (One sees in Figs. 1 and 2 that this 
factor is about 2 for the point split charge operator.) I am assuming that the 
reduction in the variance is approximately linear in subtraction order S. S 
is imagined to be a continuously variable quantity. The total time per gauge 
field is 



Tjv = {MTm + STs), (35) 

which is kept fixed as M^ff is varied. Tm is the noise time overhead and Ts 
is the subtraction time overhead. The optimum choices for S and M are now 

resulting in 

AsTl TmX X'^Ts , , 

Meff = ^ + — ^. (38) 

ATsTm 2Tm ^AsTm ^ ' 

The interesting aspect of this last result is that the effective number of noises, 
M^ff, is now quadratic in T^. This is a consequence of our observation that 
the slopes in Figs. 1 and 2 are approximately linear in S. The immediate 
implication is that for large Tjv the noise error bar can be made to vanish 
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like the inverse of the simulation time rather than as the usual inverse square 
root, at least in the range of the existing linear behavior. 

These equations are also helpful when one has an exact analytical rep- 
resentation of the trace of the perturbative series up to some order, S exact, 
making Ts zero. Then, by comparing the Me// values in the two cases, one 
may find the lowest value of S, Slowest, such that M^jy* > M^ff*. For 
example, when X f« and for common values of Tjv and Tm in the two 
simulations, one has 

Slowest — "^Sexact- ("^9) 

That is, for the extra Monte Carlo overhead to pay off, one must attempt to 
subtract to at least twice as high an order in k as the exact evaluation. Since 
subtraction is always exact for nonlocal operators up to second order and for 
local ones up to third order, making Slowest > 6 will usually result in a more 
efficient simulation than the default exact one. 



6 Summary 

We have seen that significant savings in computer resources may be obtained 
by partitioning the Wilson matrix appropriately. An efficient partitioning 
reduces the variance of an operator by leaving out the largest off-diagonal 
matrix elements of the quark propagator so that in spite of having to do more 
inversions, a smaller variance is produced for the same amount of computer 
time. A numerical investigation in Dirac and color spaces revealed efficient 
partitionings for 24 local and nonlocal operators summarized near the end of 
Section 3. 

We have also seen that large time savings are possible using subtraction 
methods for selected operators in the context of unpartitioned noise simula- 
tions. This method was shown to be effective for the scalar and local vector 
currents, but most effective for the point-split vector currents. Since pertur- 
bative subtraction is based on the hopping parameter expansion of the quark 
propagator, such methods can become less effective at lower quark masses, 
although we found the variance reduction was still quite significant for the 
point-split vector operator at k = 0.152. Similar methods can be devised for 
other operators (axial, pseudoscalar, tensor) by implementing these ideas in 
the context of Dirac/color partitioned noise methods. 

There are still a number of open questions here. The reasons for the 
strange partitioning patterns found in Section 3 are not known. In addi- 
tion, the reason why the variance of some operators respond more sensitively 
to perturbative subtraction than others is obscure. These questions are im- 
portant because their answers could lead one to better simulation methods. 
Another question is how far the linear slopes in Figs. 1 and 2 persist at high 
subtraction orders. Since the perturbative expansion of the Wilson matrix 
does not converge at small quark mass, the slope of such curves probably 
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levels off at high enough S. We have seen, however, that before this leveling 
off occurs the number of effective noises grows quadratically in the simula- 
tion time. It was pointed out that this implies that the noise error bar can 
be made to vanish like the inverse of the simulation time in the range of the 
existing linear behavior. 
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